10 - Windpark Essertine-Sur-Rolle

In diesem Skript wird eine Abschätzung des Windpotentials für Rolle und Mont-sur-Rolle gemacht. Es ist bekannt, dass in Essertine-Sur-Rolle ein Windpark für den Jahresbedarf von ca. 5400 Haushalte entstehen soll. Folgendes Vorgehen wird verfolgt:

In einem ersten Schritt prüfen wir, ob Rolle und Mont-Sur-Rolle im Perimeter der nächsten 5400 Haushalte zu liegen kommen.
Im folgenden wird angenommen, dass die beiden Gemeinden einen Anteil der gesamten Produktion des Windkraftwerkes zur Verfügung haben, welcher gerade dem Anteil an den insgesamt 5400 Haushalten entspricht (n Haushalte Roll / 5400).
Nun wird ein gemessenes Windprofil (Stundenwerte 2017) der Station Saint-Prex (nächste Windstation von Meteo Schweiz) so aufbereitet, dass nur für die Stromproduktion relevante Geschwindigkeiten vorkommen. Zudem wird das Windprofil geglättet.
Das Jahrepotential Wind, welches Rolle und Mont-Sur-Roll vom neuen Windpark zusteht, wird nun proportional auf das Windprofil aufgeteilt, so dass nun Stundenwerte zur verfügung stehen und die Summe noch immer dem Jahrepotential entspricht.

In einem spätere Schritt kann dieses Potential-Profil, entsprechend der Bedarfsprofile, auf jeden einzelnen Haushalt übertragen werden.

Daten einlesen

Die Daten aus dem vorhergehenden Arbeitsschritt werden eingelesen. Hinzu kommt das Windprofil von Saint-Prex und die GWS Daten.

load("output/05_output/geom/05_rolle_bld_out.Rdata")
rolleCRS<-proj4string(rolle_bld_out)

#gws von Rolle und Umgebung einlesen
rolle_gws<-readOGR(layer = "rolleUmgebung", dsn = paste(dataPath,"18_GWS/",sep=""))

## OGR data source with driver: ESRI Shapefile 
## Source: "/Volumes/data$/ta/60 FuE/6096 SCCER/609635 FEEBD-II/60963505 Work Packages/JA RED/06-Daten/18_GWS", layer: "rolleUmgebung"
## with 4667 features
## It has 161 fields
## Integer64 fields read as strings:  G13TOT G13A01 G13A02 G13A03 G13A05 G13A06 G13A07 G13A08 G13B01 G13B02 G13B03 G13B04 G13B05 G13B06 G13B07 G13B08 G13B09 G13B10 G13B11 G13G01 G13G02 G13G03 G13G04 G13G05 G13G06 G13G07 G13G08 G13G09 G13H01 G13H02 G13H03 G13H04 G13H05 G13H17 G13H19 G13E01 G13E02 G13E03 G13E04 G13E05 G13E06 G13E07 G13E08 G13E17 G13E19 G13EWW01 G13EWW02 G13EWW03 G13EWW04 G13EWW05 G13EWW06 G13EWW07 G13EWW08 G13EWW17 G13EWW19 G13W00 G13W01 G13W02 G13W03 G13W04 G13W05 G13W06 G13W07 G13W10 W13TOT W13TF W13T01 W13T02 W13T03 W13T04 W13T05 W13T06 W13T07 W13T08 W13T09 W13T10 W13T11 W13T1 W13T1F W13T101 W13T102 W13T103 W13T104 W13T105 W13T106 W13T107 W13T108 W13T109 W13T110 W13T111 W13T2 W13T2F W13T201 W13T202 W13T203 W13T204 W13T205 W13T206 W13T207 W13T208 W13T209 W13T210 W13T211 W13T3 W13T3F W13T301 W13T302 W13T303 W13T304 W13T305 W13T306 W13T307 W13T308 W13T309 W13T310 W13T311 W13T4 W13T4F W13T401 W13T402 W13T403 W13T404 W13T405 W13T406 W13T407 W13T408 W13T409 W13T410 W13T411 W13T5 W13T5F W13T501 W13T502 W13T503 W13T504 W13T505 W13T506 W13T507 W13T508 W13T509 W13T510 W13T511 W13T6 W13T6F W13T601 W13T602 W13T603 W13T604 W13T605 W13T606 W13T607 W13T608 W13T609 W13T610 W13T611

rolle_gws<-spTransform(rolle_gws,CRSobj = rolleCRS)

rolle_gws@data<-rolle_gws@data[,c("OBJECTID","W13TOT","W13TF")]

#ort des windkraftwerks
wind_loc<-SpatialPointsDataFrame(coords = matrix(data = c(512402.22, 148480.89),nrow = 1), 
                                 data = data.frame(id=1),
                                 proj4string = CRS(rolleCRS))

#windprofil st-prex
wind<-read.csv(paste(dataPath,"14_wind/stPrex_10minMean_wind_2017.txt",sep=""),sep=";", na.strings = "-")
names(wind)<-c("st","time","wind")

wind$wind<-wind$wind/3.6

Abdeckung Windpark

Die Abdeckung des Windparks wird anhand des Wissens um ein Potential von 5400 Haushalten und mit den Daten des GWS geschätzt. Die GWS Daten enthalten als Attribut die Anzahl Wohnungen pro Hektar.

Verschneidung der GWS Daten mit dem Einflussgebiet des Windparks

In einem ersten Schritt wird die Distanz vom geplanten Standort des Windparks zum Zentrum jedes Hektarrasters des GWS berechnet. In Folge werde die Anzahl in Reihenfolge zunehmender Distanz summiert.

gws_dist<-spDistsN1(coordinates(rolle_gws),wind_loc,longlat = F)
rolle_gws$distToWind<-gws_dist

rolle_gws_dat<-rolle_gws@data
rolle_gws_dat<-rolle_gws_dat[order(rolle_gws_dat$distToWind,decreasing = F),]

rolle_gws_dat$W13TOT<-as.numeric(as.character(rolle_gws_dat$W13TOT))
rolle_gws_dat$W13TF<-as.numeric(as.character(rolle_gws_dat$W13TF))

rolle_gws_dat$W13TOT_distSum<-rolle_gws_dat$W13TOT
rolle_gws_dat$W13TF_distSum<-rolle_gws_dat$W13TF
for(i in 2:nrow(rolle_gws_dat)){
  rolle_gws_dat$W13TOT_distSum[i]<-rolle_gws_dat$W13TOT[i]+rolle_gws_dat$W13TOT_distSum[i-1]
  rolle_gws_dat$W13TF_distSum[i]<-rolle_gws_dat$W13TF[i]+rolle_gws_dat$W13TF_distSum[i-1]
}

rolle_gws_dat_orig<-rolle_gws@data
rolle_gws_dat_orig<-rolle_gws_dat_orig%>%
  left_join(rolle_gws_dat[,c("OBJECTID","W13TOT_distSum","W13TF_distSum")],by="OBJECTID")

rolle_gws@data<-rolle_gws_dat_orig

Visualisierung

In folgender Karte werden einerseits die Summen von Wohnungen (in Abhängigkeit der Distanz zum Windpark) gezeigt. Zudem wird das maximale Potential von 5400 Wohnungen räumlich dargestellt.

wind_loc.wgs<-spTransform(wind_loc,CRSobj = CRS("+init=epsg:4326"))
rolle_gws.wgs<-spTransform(rolle_gws,CRSobj = CRS("+init=epsg:4326"))
rolle_bld_out.wgs<-spTransform(rolle_bld_out,CRSobj = CRS("+init=epsg:4326"))

pal <- colorBin(
  palette = "YlGnBu",
  domain = rolle_gws.wgs$W13TOT_distSum,
  bins=quantile(rolle_gws.wgs$W13TOT_distSum,na.rm=T,seq(0,1,length.out = 8))
)

qpal_bed <- colorBin(
  palette = "Reds",
  domain = rolle_bld_out.wgs$enerBed_awel,
  bins=quantile(rolle_bld_out.wgs$enerBed_awel,na.rm=T,seq(0,1,length.out = 8))
)

#visualisierung mit leaflet
m <- leaflet() %>%
    
  addProviderTiles(providers$OpenStreetMap, group = "normal") %>%
  addProviderTiles(providers$Stamen.TonerLite, group = "OSM (b/w)") %>%
  addProviderTiles("Esri.WorldImagery", group = "ESRI") %>%
  
  addPolygons(data=rolle_gws.wgs,
              stroke = TRUE,
              fillOpacity = 0.8, 
              color = ~pal(rolle_gws.wgs$W13TOT_distSum),
              weight=1,
              group = "Anzahl Wohnungen (GWS)")%>%
  
  addPolygons(data=rolle_gws.wgs[rolle_gws.wgs$W13TOT_distSum<5400,],
              stroke = TRUE,
              fillOpacity = 0.5, 
              color = "grey30",
              weight=1,
              group = "Haushalte im Umkreis (< 5400)")%>%
  
  addCircleMarkers(lng = wind_loc.wgs@coords[,1],
                   lat = wind_loc.wgs@coords[,2],
              fillOpacity = 0.8, 
              color = "red",
              weight=1,
              group = "Standort Windkraftwerk")%>%
  
  addPolygons(data=rolle_bld_out.wgs,
              stroke = TRUE,
              fillOpacity = 0.8,
              color = ~qpal_bed(rolle_bld_out.wgs$enerBed_awel),
              weight=1,
              group = "Energiebedarf (AWEL)")%>%
  
  addLayersControl(
    baseGroups = c("OSM (b/w)","normal","ESRI"),
    overlayGroups = c("Anzahl Wohnungen (GWS)","Haushalte im Umkreis (< 5400)","Standort Windkraftwerk","Energiebedarf (AWEL)"),
    options = layersControlOptions(collapsed = F)
  )%>%
  addLegend(pal = pal, values = rolle_gws.wgs$W13TOT_distSum,
    title = "Anzahl Wohnungen")%>%
  addLegend(pal = qpal_bed, 
            values =rolle_bld_out.wgs$enerBed_awel, 
            opacity = 1,
            title = "Energiebedarf: kWh/a",
            position = "bottomleft"
            )

#leaflet karte ausführen
m

Die Karte zeigt, dass die meisten Häuser in Roll Teil der 5400 Haushalte sind, welche am nächsten zum neuen Windpark in Essertine-Sur-Rolle liegen. Mont-Sur-Rolle ist ganz im Perimeter enthalten. Es ist daher anzunehmen, dass ganz Rolle erschlossen werden wird.

Windpark-Potential als Stundenprofil

Das Jahrespotential des Windparkes wird auf die Haushalte in Rolle heruntergerechnet. Das entprechende Potential wird in ein Jahresprofil mit Stundenwerten umgerechnet. Diese Umrechnung wird proportional entlang eines Windprofils einer nahe gelegenen Windstation gemacht (Saint-Prex).

Explorative Windanalyse

In einem ersten Schritt wird das Windprofil von Saint-Prex dargestellt.

year<-substr(as.character(format(wind$time, scientific=F)),start = 1,stop = 4)
month<-substr(as.character(format(wind$time, scientific=F)),start = 5,stop = 6)
day<-substr(as.character(format(wind$time, scientific=F)),start = 7,stop = 8)
hour<-substr(as.character(format(wind$time, scientific=F)),start = 9,stop = 10)
min<-substr(as.character(format(wind$time, scientific=F)),start = 11,stop = 12)
date<-paste(year,".",month,".",day," ",hour,":",min,sep="")
wind$date<-as.POSIXct(date, format="%Y.%m.%d %H:%M", tz = "GMT")

wind.agg<-wind%>%
  mutate(index = as.POSIXct(round(date,"hours")),
         wind=as.numeric(as.character(wind)),
         date=NULL,
         time=NULL,
         st=NULL)%>%
  dplyr::group_by(index) %>%
  dplyr::summarise(windH=mean(wind,na.rm=T)) %>%
  mutate(hour=as.numeric(strftime(index, format="%H")),
                weekDay=lubridate::wday(as.Date(index),week_start=1),
                week = as.numeric(format(as.Date(index),"%W")),
                month = as.numeric(format(as.Date(index),"%m")),
                season = cut(month,breaks = seq(0,12,by = 2),labels = c(1,2,3,3,4,1)),
                year = as.numeric(format(as.Date(index),"%Y")),
                weekYear = paste(week,year,sep="_"),
                monthYear = paste(month,year,sep="_"))

levels(wind.agg$season)<-c("Winter","Spring","Summer","Autumn")

ggplot(wind.agg)+
  geom_path(aes(x=index,y=windH),color="black")+
  theme_minimal()+
  ggtitle("Windprofil 2017")+
  ylab("Windgeschwindigkeit (m/s)")+
  xlab("Tageszeit")

Das Jahresprofil wird im folgenden noch auf Saisonalitäten untersucht.

p3<-wind.agg%>%
  dplyr::group_by(season,hour)%>%
  dplyr::summarise(windMean=mean(windH,na.rm=T),
            windSD=sd(windH,na.rm=T))%>%
  ggplot(.)+
  geom_path(aes(x=hour+1,y=windMean),color="black")+
  geom_path(aes(x=hour+1,y=windMean-windSD),color="grey80")+
  geom_path(aes(x=hour+1,y=windMean+windSD),color="grey80")+
  facet_wrap(.~season)+
  theme_minimal()+
  ggtitle("Tagesverlauf pro Jahreszeit")+
  ylab("Windgeschwindigkeit (m/s)")+
  xlab("Tageszeit")+
  scale_x_continuous(breaks = 1:24)

p3

Die vier Jahreszeiten zeigen stark unterschiedlichen Tagesverläufe. Am Nachmittag ist jeweils ein leichter Anstieg des Windes zu verzeichnen. Die Windgeschwindigkeit ist mit ca. 3-4 m/s über das ganze Jahr hinweg relativ konstant und stark.

Windpotential für Rolle

Die oben besprochene Umrechnung wird vorgenommen. Das Jahrespotential der Windanalage wird auf den Anteil der Häuser in Rolle heruntergebrochen. Das Potential, welches Rolle zur Verfügung steht, wird auf das Windprofil übertragen. Die Jahressumme bleibt dabei erhalten, wird aber proportional auf Stundenwerte übertragen.
Für diese Berechnung wird das Windprofil erst leicht angepasst. Zum einen werden die Spitzen-Windgeschwindigkeiten von über 10 m/s auf konstant 10 m/s heruntergesetzt. Werte tiefer als 2 werden auf 0 gesetzt, da bei solch leichtem Wind keine Elektrizität produziert werden kann. Die Werte sind einer typischen Kennlinie von Windkraftanalgen entnommen. Als letzte Anpassung werden die Stundenwerte leicht geglättet, damit die Werte ein bisschen weniger Rauschen und die Partikularitäten des Jahres 2017 weniger zu Tage treten.

#anzahl haushalte in Rolle aus GWS
rolle.ch<-gConvexHull(rolle_bld_out)
ov<-over(rolle_gws,rolle.ch)
rolle_gws.s<-rolle_gws[!is.na(ov),]
nHaus<-sum(as.numeric(as.character(rolle_gws.s$W13TOT)))

#potential für Rolle als Anteil der Rolle Haushalte an allen 5400 Haushalte, die beliefert werden
potRolle<-nHaus/5400*(5400*3500)

#Windprofil einschränken
#<2 -> 0
#>15 -> 15
wind.agg<-wind.agg%>%
  mutate(windEing = if_else(windH<2, 0,if_else(windH>=10, 10,windH)),
         windEing.ma = rollmean(x=windEing, 5, align="center",fill = NA))

wind.agg$windEing.ma[is.na(wind.agg$windEing.ma)]<-mean(wind.agg$windEing.ma,na.rm=T)

wind.sum<-sum(wind.agg$windEing,na.rm=T)

potRolle.rel<-potRolle/wind.sum

wind.agg$potRolleRel<-potRolle.rel*wind.agg$windEing

ggplot(wind.agg)+
  geom_path(aes(x=index,y=potRolleRel),color="black")+
  theme_minimal()+
  ggtitle("Jahresprofil Windpotential")+
  ylab("Windpotential Rolle (kWh)")+
  xlab("Tageszeit")

Das resultierende Profil zeigt das Windpotential, welches Rolle über den Verlauf eines Jahres (pro Stunde) zur Verfügung steht. Wenn beispielsweise die Bedarfsprofile für Elektrizität für alle Gebäude von Rolle bekannt sind (bzw. geschätzt wurden), kann das Windpotential proportional auf die Gebäude aufgeteilt werden. Damit ist eine Schätzung gegeben, wie viel des Elektrizitätsbedarf pro Stunde und Gebäude durch Wind abgedeckt werden kann.

Resultate speichern

Alle wichtigen Informationen werden abschliessend gespeichert…

#html karte speichern
wd<-getwd()
htmlwidgets::saveWidget(m, 
                        file=paste(wd,"/output/10_output/map/10_windMap.html",sep=""),
                        selfcontained = T)

#windprofil speichern
save(wind.agg,file="output/10_output/geom/10_windProfile.Rdata")